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Abstract 

In this work is obtained an explicit form for the perturbation effect on the matrix of regression 
coefficients on the optimal solution in multiresponse surface methodology. Then, the sensitivity 
analysis of the optimal solution is studied and the critical point characterisation of the convex 
program, associated with the optimum of a multiresponse surface, is also analysed. Finally, the 
asymptotic normality of the optimal solution is derived by standard methods. 
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1 Introduction 

The multiresponse surface methodology has been considered as a very useful tool in the study 
of designs, phenomena and experiments. Which enables us to propose a set of analytical 
relationship between responses and controlled variables through a process of continuous 
improvement and optimisation. 

It is assumed that a researcher knows a system and a corresponding set of observable 
responses variables Y\, ■ ■ ■ ,Y r which depends on some input variables, xi, . . .x n . This also 
supposes that that the input variables x^ can be controlled by the researcher with a mini- 
mum error. 

Typically we have that 

ifc( x ) =Vk{xu---X n ), k = 1, . . . ,r, and x = (x u . . .x n )', (1.1) 
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where the form of the functions ??fc(-)'s are unknown and perhaps, very complex, and it is usu- 
ally termed as the true response surface. The success of the response surfaces methodology 
depends on the approximation of %(•) for a polynomial of low degree in some region. 

For purposes of this paper is assumed that %(•) can be soundly approximated by a 
polynomial of second order, that is 



r ft (x) = /? ofc + ^A 



kXi 



i=l 



(1.2) 



i=l j>i 



where the unknown parameters /3 - s can be estimated via regression's techniques, as it will 
be described in next section. 

Next, we are interested in obtaining the levels of the input variables x? such that the 
response variables Yi, • • • ,Y r are simultaneously minimal (optimal). This can be achieved 
if the following multiobjetive mathematical program is solved 

/ n(x) \ 



nun 



V ^(x) / 
subject to 

xe x, 



(1.3) 



where X is certain operating region for the input variables 
Now, two questions, closely related, can be observed: 



1. 



When the estimations of (|1.2|) for k = 1, . . . ,r arc considered into (|1.3jl . the critical 
point x* obtained as solution shall be a function of the estimators s of the (3j s . 
Thus, given that /3j 3 are random variables, then x* = x*(/3j s ) is a random vector 

too. So, under the assumption that the distribution of /3 is known, then, what is the 
distribution of x*(/3.' s )? 



2. And, perhaps it is not sufficient to know only a point estimate of x*(/3 J s ), could be 
more convenient to know a estimated region or a estimated interval. 

In particular, t he distribution of the critical poi n t in a univa riate response surface model 
was studied bv iDiaz Garcia and Ramos- Quiroera ( 2001 . 2002 ). when y(x) is defined as an 
hypcrplane. 

Now, in the context of the mathematical programming problems, the sensitivity analysis 
studies the effect of small perturbations in: (1) the parameters on the optimal objective 
function value and (2) the critical point. In general, these parameters shape the objective 
func tion and const r aint t h e approach to the mat hematical programming pr oblem. In partic- 
ular, Jagannathanl ( 19771 ). Dupacova (1984) and Fiacco and Ghaemi ( 19821 ) have studied the 
sensitivity analysis of the mathematical programming, among many other authors. As an 
immediate consequence of the sensitivity analysis emerges the asymptotic normality study 
of the critical point, which can be performed by standard methods of mathematical statis- 
tics (s ee similar results for the case of maximum likelihood estimates lAitchison and Silvevl 
(1958)). This last consequence makes the sensitivity analysis very appealing for researching 
from a statistical point of view. However, this approach must be fitted into the classical 
philosophy of the sensitivity analysis; i.e., we need to translate into the statistical lan- 
guage, the general sensitivity analysis methodology, which deals with a number of ways in 



Asymptotic normality in multirespon.se surface 



3 



which the estimators of certain model are affected by omission of a particular set of vari- 
ables or by the inclu s ion or omission of a particular observation or set of observations, see 
Chatteriee and Hadil (jl988l) . 

This papers pursues to important aims: the effect of perturbations of the matrix of 
regression parameters on the optimal solution of the multiresponse surface model and the 
asymptotic normality of the critical point. First, in Scction[3]somc notation is defined. Then, 
the multiresponse surface mathematical program is proposed in Scction|3]as a multiobjective 
mathematical programming problem and a general solution is considered in terms of a 
functional. The characterisation of the critical point is given in Section |4] by stating the 
first-order and second-order Kuhn- Tucker conditions. Finally, the asymptotic normality of 
a critical point is established in Section [5] and for a particular form of the functional, the 
asymptotic normality of a critical point is also derived. 



2 Notation 

For convenience, the principal properties and usual notations are give n here. A detailed 
discus sion of the m ultiresponse surface meth odology can be found in iKhuri and Cornell 
(Il987l . Chap. 7) and lKhuri and Conlonl (Il98ll) . 

Let N be the number of experimental runs and r be the number of response variables, 
which can be measured for each setting of a group of n coded variables xi, 22, • ■ • , x n . It is 
assumed that the response variables can be modelled by a second order polynomial regression 
model in terms of Xi, i = 1, . . . , n. Hence, the k th response model can be written as 



Xfc/3, 



1. 



. , r, 



(2.4) 



where Y& is an N x 1 vector of observations on the k th response, X& is an N x p matrix of 
rank p termed the design or regression matrix, p = 1 + n + n(n + l)/2, (3 k is a p x 1 vector 
of unknown constant parameters, and is a random error vector associated with the k th 
response. For purposes of this study is assumed that Xi = • • • = X r = X. Therefore, (|2.4[) 
can be written as 

Y = XB + E (2.5) 



where Y 



y i: y 2 



:Y r 



/3ii/3 2 : • • • & 



, moreover 



Pk = (Pok,Pik, ■ ■ -,Pnk,P 



life, 



, Pnnki Pl2k, ■ ■ ■ 1 P(n-l)nk) 



and E 



ei:e 2 : 



:e r , such that E ~ J\f Nxr (0,I N <g> S) i.e. E has an N x r matrix 
multivariate normal distribution with E(E) = and Cov(vecE') = Ijy (g) S, where S is a 



r x r positive definite matrix. Now, if A 



Ai:A 2 



with A,-, j 



1, • ■ ■ , r the 



columns of A; then vec A = (A\ , A9, . . . , AD' and (g> denotes the direct (or Kronecker) 
product of matrices, see iMuirheadl j 19821 Theorem 3.2.2, p. 79). In addition denote 



- x = [xi, X2, ■ ■ ■ , x n ) : The vector of controllable variables or factors. Formally, an Xi 
variable is associated with each factor A, B, ... 



-- / 9 1 :/3 2 :---:/3 
where 



The least squares estimator of B given by B = (X'X) 'X'Y, from 
Pk = (X'X) _1 X'Y fe = (/3 fe, /3ifc, • • • , PmPllk, ■ ■ ■ ,0nnk,Pl2k, P{n-l)nkY 
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k = 1,2, ... , r. Moreover, under the assumption that E ~ A/jvxr(0, Ijv ® E), then 
B - Af pxr (®, (X'X)" 1 ® £), with Cov(veci') = (X'X)" 1 ® S. 

- z(x) = (1, Xi, X2, • ■ • 5 X n , X-^, • ■ • 7 X\X2i X\X% . . . , X n _iX n ) . 

- 3ift = (Pik,---,Pnk)' and 

^ 2/3llfc /3i2fc • • • Pink ^ 
/3aifc 2/3 2 2fc ■ • ' Pink 



B k = 



Y fc (x) 



A, 



2 k 



^Pnnk J 



pok+^2PikXi+^2p 



iik^i ~t~ ^ ^ ^ ^ Pijk%i-Ej 
i—1 i—1 j>i 



i=l 

= p ok + 3ifeX + x B fe x : 
The response surface or predictor equation at the point x for the k" 1 response variable. 



- Y(x) = ^Fi(x), F 2 (x), . . . , F r (x)J = B'z(x): The multiresponse surface or predicted 

response vector at the point x. 

- Y'(Ijv - X(X'X)- 1 X')Y , r , . ^ 

- Yj = : The estimator of the variancc-covariancc matrix £ such 

N —p 

that (TV — p)H has a Wishart distribution with (N — p) degrees of freedom and the 
parameter S; this fact is denoted as (N — p)S ~ W r (N — p, £). Here, I m denotes an 
identity matrix of order m. 



- Finally, note that 
and 



E(Y(x)) = £(B'z(x)) = B'z(x) 



Cov(Y(x)) = z'(x)(X'X)- 1 z(x)S. 
An unbiased estimator of Cov(Y(x)) is given by 

Cov(Y(x)) = z'(x)(X'X)- 1 z(x)S. 



(2.6) 
(2.7) 

(2.8) 



3 Multiresponse surface mathematical programming 

In the following sections, we make use of the multiresponse mathematical programming and 
multiobjective mathematical programming. For convenience, the concepts and notations 
required are listed below in terms of the estimated model of multirespo nse surface mathe- 
matic al programming . Definition s and detailed prop e rties can be fo und in Khuri and Conlon 
(|l98lh . lKhuri and Cornell (|l987t) . iRfos et all (fl989l) . ISteuerl (jl986l) . and iMiettinenl (|l999f ). 

The multiresponse mathematical programming or multiresponse optimisation (MRO) 
problem is proposed, in general, as follows 



min Y(x) = min 



/ n(x) \ 

Y 2 (x) 



V Y r (x) / 



(3.9) 



subject to 
x £ X. 
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It is a nonlinear mult iobjective mathe matical programming problem, see Steuer ( 19861 ). 
Rios et all (|l989h and iMiettinenl (|l999f ); and X denotes the experimental region, usually 
taken as a hypercube 

X = {x\k < Xi < Ui, i = 1, 2, . . . , n}, 

where 1 — (l\, I2, ■ ■ . , l n )' , defines the vector of lower bounds of factors and u = {u\, U2, ■ ■ ■ , u n ) 
gives the vector of upper bounds of factors. Alternatively, the experimental region can taken 
as a hypersphere 

£ = {x|x'x< c 2 ,c£ 5ft}, 

where , c is set according to the experimental design model under consideration, see lKhuri and Cornelll 
(|l987t ). Alternatively fl3U) can be written as 



min Y(x). 

In the response surface methodology context, the multiobjective mathematical programs 
rarely contains a point x* which can be considered as an optimum, i.e. few cases satisfy 
the requirement that Yfc(x) is minimum for all k = 1,2, ... ,r. From the viewpoint of 
multiobjective mathematical programming, this justifies the following notion of the Pareto 
point: 

We say that Y*(x) is a Pareto point of Y(x), if there is no other point 
Y x (x) such that Y x (x) < Y*(x), i.e. for all k, Y£(x) < Y*(x) and 
Y*(x) ^ Y*(x). 

Steuer ( 19861 ). Illios et all ( 19891 ) and Miettinen ( 19991 ) established the existence criteria for 
Pareto points in a multiobjective mathematical programming problem and the extension of 
scalar mathematical programming (Kuhn-Tucker's conditions) to the vectorial case. 

Methods for solving a multiobjective mathematical program are based on the existing 
information about a particular problem. There are three possible sc enarios: w h en th e 
investigator poss esses either compl ete, partial or null information, see iRios et al. (|l989l) . 
Miettinenl (|l999t ) and ISteuerl (|l986l ). In a response surface methodology context, complete 
information means that the investigator understands the population in such a way that it is 
possible to propose a value function reflecting the importance of each response variable. In 
partial information, the investigator knows the main response variable of the study very well 
and this is sufficient support for the research. Finally, under null information, the researcher 
only possesses information about the estimators of the response surface parameter, and with 
this elements an appropriate solution can be found too. 

In general, an approach for solving a multiobjective mathematical program consist of 
proposing an equivalent nonlinear scalar mathematical program, i.e. as a solution of 
is proposed the following problem 



min / ^Y(x) 

subject to 
x £ X, 



(3.10) 



where /(•) defines a functional (/(•) is a function that takes functions as its argument, i.e. 
a function whose domain is a set of functions). Moreover, in the context of multiobjective 
mathematical programmming, the functional /(•) is such that if 9Jt C W denotes a set of 
multiresponse surface functions, then 

The functional is a function f : 9H — > 5ft such that minY(x*) < 
minY( Xl ) ^ /(Y(x*)) < /(Y( Xl )), x'/xj. 
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In order to consider a greater number of potential solutions of (|3.9p . usually studied in the 
multicriteria mathematical programming, the following alternative problem to p,10[) can be 
proposed 



/(Y(x) 



subject to (3-11) 
x e xn 6, 

where & is a subset generated by additional potential constraints, generally derived by a 
particular technique used for establishing the equivalent scalar mathematical program (|3.10j) . 
In some particular cases of p.!0[) . a new fixed parameter may appear, a vector of response 
weights w = (wi, u>2, ■ ■ ■ , w r )' , and/or a vector of target values for the response vector r = 
( T i: T 2, ■ ■ ■ j T r ) . Particular examples of this equiv alent univariate objective mathematical 



programming are the use of goal programming, see iKazemzadeh et al\ (|2008l ). and of the e- 
constraint model, see iBilesI f|l975h . among many others. In particular, under the e-constraint 
model, (|3.11j) is proposed as 

min Yj (x) 

x J 

subject to 

?i(x) <n 

g-xCx)^.! (3-12) 

Y r (x) < T r 
x £ X. 

4 Characterisation of the critical point 

In the rest of the paper we shall develop the theory of the problem (|3.10[) ; it is easy to see 
that this problem can be extended with minor modifications to the problem p. lip . 

Let x*(B) € 5R™ be the unique optimal solution of program (|3.10[) with the corresponding 
Lagrange multiplier A* (B) £ 5ft. The Lagrangian is defined by 

L(x,A;B) = /(y(x)) +A(||x|| 2 -c 2 ). (4.13) 

Similarly, x* (B) £ 5ft" denotes the unique optimal solution of program (|1.3[) with the corre- 
sponding Lagrange multiplier A* (B) £ 5ft. 

Now we establish the local Kuhn- Tucker conditions that guarantee that the Kuhn- Tucker 

^ i- ^ i ' 

point r* (B) = x* (B), A* (B) £ 5ft n+1 is a unique global minimum of convex program (|3.10j) . 

df 

First recall that for / : 5ft™ — > 5ft, — = V x denotes the gradient of function /. 

ax 

Theorem 4.1. The necessary and sufficient conditions that a point x*(B) £ 5ft™ for arbitrary 
fixedM £ 5ft p , be a unique global minimum of the convex program IS.lty) is that, x*(B) and the 
corresponding Lagrange multiplier A* (B) £ 5ft, fulfill the Kuhn- Tucker first order conditions 



V x L(x,A;B) = V x /(Y(x)J + 2A(B)x = (4.14) 

V A L(x,A;B) = ||x|| 2 - c 2 < (4.15) 

A(B)(||x|| 2 -c 2 ) = (4.16) 

A(B) > (4.17) 



Asymptotic normality in multirespon.se surface 



7 



In addition, assume that strict complementarity slackness holds at x*(B) with respect 
to A*(B), that is 

A*(B) >0^||x|| 2 -c 2 =0. (4.18) 
Analogously, the Khun- Tucker condition (|4.14[) to (|4.17l) for B = B are stated next. 

Corollary 4.2. The necessary and sufficient conditions that a point x*(B) G 5ft™ for arbi- 
trary fixed B G K p , be a unique global minimum of the convex program (1.3]) is that, x*(B) 
and the corresponding Lagrange multiplier A*(B) g 5R, fulfill the Kuhn- Tucker first order 
conditions 



V x L(x,A;B) = V x /(Y(x)J + 2A(B)x = (4.19) 

V A L(x,A;B) = ||x|| 2 - c 2 < (4.20) 

A(B)(||x|| 2 -c 2 ) = (4.21) 

A(B) > (4.22) 

and A(B) = when ||x|| 2 - c 2 < at [x*(B), A*(B)]'. 

Observe that, due to the strict convexity of the constraint and objective function, the 
second-order sufficient condition is evidently fulfilled for the convex program (|3.10p . 

The next result states the existenc e of a once continuously differentiable solution to 
program (|3.10|) . see lFiacco and Ghaemil (|l982l) . 

Theorem 4.3. Assume that 18$ holds and the second-order sufficient condition is satis- 
fied by the convex program (3.10]) . Then 

1. x*(B) is a unique global minimum of program (1.3$ and A*(B) is also unique. 

2. ForM G 14(B) (is an e— neighborhood or open ball), there exist a unique once contin- 
uously differentiable vector function 



r*(B) 



x*(B) 



G 5ft 



n+l 



satisfying the second order sufficient conditions of problem (1.3]) . such that r*(B) = 
[x*(B), A*(B)]' and hence, x*(B) is a unique global minimum of problem (3.10\) with 
associated unique Lagrange multiplier A*(B). 

3. ForM G 14(B), the status of the constraint is unchanged and A*(B) > ||x|| 2 -c 2 = 
holds. 



5 Asymptotic normality of the critical point 

This section considers the statistical and mathematical programming aspects of the sensi- 
tivity analysis of the optimum of a estimated multiresponse surface model. 

Theorem 5.1. Assume: 

1. For any B G V e (M), the second-order sufficient condition is fulfilled for the convex 
program (3.10)) such that the second order derivatives 

<9 2 £(x,A;B) <9 2 L(x,A;B) <9 2 L(x,A;B) d 2 L(x,A;B) a 2 L(x,A;B) 
<9x<9x' ' dxdvec'B ' dxdX ' dXdx' ' dXdvccM 
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exist and are continuous in 



is positive definite. 



x*(B),A*(B) e T4([x*(B), A*(B)]') and 

<9 2 £(x,A;B) 
9x<9x' 



2. M v , the estimator of the true parameter vector B„, is based on a sample of size N v 
such that 

- B„) ~ AA pxr (B, 0), -^0 = (X'X)- 1 <g> S. 



5. f£7gp ZioZds /< 



or 



TTien asymptotically 



x*(B)-x*(B) 4Af„(0„,3), 



where the n x n variance- covariance matrix 



<9x* 



dvecl 



& 



<9x*(JB 
9 vec] 



As 



■0 = (X'X) -1 ® s 



such that all elements of ( <9x*(B)/<9vecB ) are continuous on any B £ 14(B); further- 
more 



<9x* 



9 vec] 



where 



[I P" 1 Q(Q'P" 1 Q)" 1 Q'] P J G, 

<9 2 L(x,A;I) 



P 

Q 



9x9x' 
9 2 L(x,A;B) 
SA<9x 



G _ d 2 L(x,\;M) 
9x9 vec' B 

Proof. According to Theorem 14.11 and Corollary 14.21 the Kuhn- Tucker conditions (14.141) 

(|4~TTj) at x*(B), A* (I) ' and the conditions (|4"H?)) - (|4T2"2"j) at [x*(B), A*(B)]' arc fulfilled for 

mathematical programs ()1.3[) and p.lO|) . respectively. From conditions ()4.19[l — f|4.22j) of 
Corollary 14. 2 1 the following system equation 



V x Z(x,A;B) = V x / Y(x 



V A L(x,A;B) = ||x|| 2 - c 2 





0, 



(5.23) 
(5.24) 



has a solution x*(B), A*(B) > 0, B. 

The nonsingular Jacobian matrix of the continuously diffcrcntiable functions (|5.23[) and 



(I5.27|) with respect to x and A at 



,A*(1 



is 



( d 2 £(x,A;B) d 2 £(x,A;B) \ 

<9x9x' dXdx 
d 2 £(x,A;B) 

\ <9x'dA 
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According to the implicit functions theorem, there is a neighborhood T4(B) such that for 
arbitrary B e V e (M), the system (|5~23| and (fOTj) has a unique solution x*(S), A*(B), I and 
by Theorem 14.31 the components of x*(B), A*(B) are continuously differcntiablc function of 
B, see 



Bigelow and Shapiro! (|1974l ). Their derivatives are given by 



/ <9x*(B) \ 



d vec! 
<9A*(1 



/ <9 2 L(x,A;B) <9 2 L(x,A;B) \ 



/ 



V d vec . 
The explicit form of (dx* 



P Q 

Q' 



dxdx' ^ 
<9 2 L(x,A;B) 



dXdx 




/ 



d 2 L(x,A ; : 

dxd vec' 1 




.(5.25) 



dx'dA 

vecB) follows from (|5 . 25[) and by the formula 



[I P 1 Q(Q'P 1 Q) 1 Q / ]P- 1 P 1 Q(Q'P 1 Q)- 1 
(Q'piQ^Q'P 1 -(Q'p-iQ)" 1 



where P is symmetric and no nsingular. 

Then from assumpti on 2, Raa (Il973l (iii), p. 388) and Bishop et al. ( 1991 . Theorem 
14.6-2, p. 493) (see also Cramerl ( 1946 . p. 353)) we have 



M n 0„, 



dx* 



dveel 



& 



/<9x*(B)Y 



vec J 



(5.26) 



Finally note that all elements of (dx* /dM) are continuous on V e (B), so that the asymptotical 
distribution (|5.26p can be substituted by 



x*(B) 



■Mn On, 



<9x*(B) 
dvccB 



& 



dx* 



d vecl 



see 



Rat] (119731 . (iv), pp.388-389). 



□ 



As a particular case, assume that the functional in (|3.10j) is defined as 

r r 

/(Y(x)) =5> fe Y fc (x), £> fc = l, 



fe=i 



fc=i 



with Wfe known constants. Then, 

Corollary 5.2. Suppose: 

1. For any B G 14(B), the second- order sufficient condition is fulfilled for the convex 
program h3.10\) such that the second order derivatives 

d 2 L(x,X:M) <9 2 L(x,A;B) <9 2 L(x,A;B) d 2 L(x,A;B) <9 2 L(x,A;B) 



dxdy 



dxd vec' 



dxdX 



dXdx' 



d\d vec 1 



exist and are continuous in 



x*(B),A* 



e K([x* 



) and 



d 2 L(x,X;M) 
dxdx' 



is positive definite. 
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2. M v , the estimator of the true parameter vector M v , is based on a sample of size N v 
such that 

1 



B„)~JV P xr(B,e), -^e = (x'x)- 1 ® s. 



3. j f^.iffp Zio/ds for B = B. T/ien asymptotically 



x*(B)-x*(B) 47V„(0„,S) 



where the n x n variance- covariance matrix 



<9vecB / \ 9vecB 



sitc/i £/ia£ elements of ^c>x*(B)/(9vecBj are continuous on any B £ 14(B); further- 
more 

V(B)x*(B)'S- 1 \ / , 
y~ V ^ - In M(x*(B) 



<9x*(B) 



dvec) 



= S" 1 



X I B s x 



where 



<9 2 L(x,A;B) 



9x<9x 



M(x) = V x z'(x 



^ =2^^B fe -2A*(B)I„. 

fe=i 

<9z'(x) 



<9x 



with 



= (0:I n :2diag(x):Cii ■ • • :C n _i) £ V""'-. 

,* = l,...,n- 1, Oj e = 1, 1; 



Oj-i 
x'Ai 

observing that when i — 1 fz.e. j — row does no£ appear in Ci; and 



A,; = 



0', 

V In-, y 



, 0' fc eH"-*,fc = !,...,». 



Proof. The required result follows from Theorem 15.11 and observing that in this particular 



case 



M(x) w kPk + 2A(B)x 



V x i(x, A;B) = < 



fc=i 



53 *"k \Pik + 2B fe x ] + 2A ( B ) X 
I, fc=i 

V Aj L(x,A;/3) = ||x|| 2 -c 2 = 



□ 
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Conclusions 



As a consequence of Theorem l4.3l now is feasib le to establish confid ence regions and intervals 



and hypothesis tests on the critical point, see lBishop et al\ (|199lL Section 14.6.4, pp. 498 



500); it is also possible to identify operating conditions as regions or intervals instead of 
isolated points. 

The results of this paper can be taken as a good first approximation to the exact problem. 
However, unfortunately in many applications the number of observations is relatively small 
and perhaps the results obtained in this work should be applied with caution. 
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